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Abstract 

For future membrane based structures, quantified 
predictions of membrane wrinkling behavior in terms of 
amplitude, angle and wavelength are needed to 
optimize the efficiency and integrity of such structures, 
as well as their associated control systems. For 
numerical analyses performed in the past, limitations on 
the accuracy of membrane distortion simulations have 
often been related to the assumptions made while using 
finite elements. Specifically, this work demonstrates 
that critical assumptions include: effects of gravity, 
supposed initial or boundary conditions, and the type of 
element used to model the membrane. In this work, a 
0.2m 1 2 membrane is treated as a structural material with 
non-negligible bending stiffness. Mixed Interpolation of 
Tensorial Components (MITC) shell elements are used 
to simulate wrinkling behavior due to a constant applied 
in-plane shear load. Membrane thickness, gravity 
effects, and initial imperfections with respect to flatness 
were varied in numerous nonlinear analysis cases. 
Significant findings include notable variations in 
wrinkle modes for thickness in the range of 50 pm to 
1000 pm, which also depend on the presence of an 
applied gravity field. However, it is revealed that 
relationships between overall strain energy density for 
cases with differing initial conditions are independent 
of assumed initial conditions. In addition, analysis 
results indicate that the relationship between amplitude 
scale (W/t) and structural scale (L/t) is linear in the 
presence of a gravity field. 

Keywords: wrinkle, simulation, finite elements, mixed 
interpolation of tensorial components (MITC), 
membrane structures 

Introduction 

With the need to develop increasingly more 
efficient and lighter mechanical and structural systems 
for spacecraft, many researchers have for decades been 
experimenting with various types of polymer 
membranes and working to simulate their behavior. 1,2 
To date, the most detailed simulations have been 
limited to focusing on global distortion and the overall 


stress field; analysis tools have not been robust enough 
to quantify wrinkle features such as amplitude, 
wavelength and orientation angle. In the future, 
quantified predictions of membrane wrinkling behavior 
in terms of these three characteristics can be used to 
optimize the efficiency and integrity of membrane 
structures and associated control systems. 

Prominent space-based examples of structural 
systems employing thin membrane materials include 
the Echo-1 passive communications satellite, high 
altitude balloons, Fresnel lens reflectors, and a sun 
shield for the James Webb Space Telescope (JWST), as 
illustrated in Figure 1 (a)-(d), respectively. Each of 
these structures is sensitive to wrinkle patterns in 
unique ways, significantly impacting the function of 
each. 



Fig. 1 Space-based membrane structures: (a) Echo- 
1, (b) high altitude balloon, (c) Fresnel lens, (d) a 
proposed sun shield concept for JVVST 


For membrane surfaces designed to be specular, 
such as Echo-1 and the lens, wrinkles cannot be 
tolerated. Wrinkles are an asset, however, to the high 
altitude balloon design. Wrinkle tolerance lying 
somewhere in between these two extremes is given by 
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the sun shield design. That is, wrinkles can be tolerated 
up until they might induce buckling in the inflated 
booms or up to the point that wrinkle amplitudes cause 
adjacent layers of membrane surfaces to come into 
contact. 

The analytical study documented herein includes 
wrinkle simulations for the basic case of a square 
membrane fixed at each comer and subjected to in- 
plane shear, which is a case that has not yet been treated 
in detail by previous works 5 . Use of the Mixed 
Interpolation of Tensorial Components (MITC) shell 
element has enabled this detailed quantitative study for 
reasons discussed in the following section. The 
subsequent section describes numerous nonlinear 
analysis cases considered, in which membrane 
thickness, gravity effects, and initial imperfections with 
respect to flatness were varied. Simulation results are 
presented and quantified in terms of wrinkle amplitude, 
orientation and wavelength. 

Results of this work are significant as a tool for: (1) 
assessing numerical analysis sensitivities in preparation 
for data correlation, and (2) developing alternative 
methods for analyzing membrane wrinkling behavior 
without relying on finite element methods. 3 

Wrinkline Theory 

Although a complete history of the development of 
theories for explaining wrinkling behavior are beyond 
the scope of this paper, three theories traditionally used 
to study and simulate wrinkled membranes are briefly 
summarized here. 

First, Tension Field Theory (TFT) mathematically 
describes wrinkling in terms of solid mechanics. 4 It 
asserts that membranes can be described as structures of 
sufficiently little thickness such that they have zero 
bending stiffness. As a result, they have no capacity to 
support compressive load or carry compressive stress. 

Plate Theory is assumed when considering 
structural responses that are dominated by in-plane 
structural stiffiiess. Out-of-plane, or transverse loads are 
considered negligible by comparison. Elements treated 
according to this theory are not expected to undergo 
significant amounts of bending, although they do have 
significant thickness. 

Classic Shell Theory is assumed when a thin 
structural element is loaded both in-plane and out-of- 
plane, and bending is expected. The theory includes 
non-negligible bending stiffness and can be used for 


5 At the time of this study, the work by X. Su, et. ai, 
entitled “Wrinkling Analysis of a Kapton Square 
Membrane Under Tensile Loading” had not yet been 
published in the proceeds from the 44* 
A1AA/ASME/ASCE/AHS Structures, Structural 
Dynamics and Materials Conference, 7-10 April, 2003. 


considering nonlinear deformation, or strain. Further, it 
considers that the shell is an inextensible surface and is 
modeled at the mid-plane of that surface - therefore, 
transverse stresses through the thickness of a surface 
are not captured. 

For elastic materials, such as membranes 
considered in this work, inextensibility is a fundamental 
principle stating that the length of a distorted path 
cannot exceed the length of an undistorted path of the 
same structure. In other words, geometric compatibility 
between initial and final states must be preserved. For 
wrinkling phenomena, this principle relates to the 
transverse direction. Hence, neither Plate nor Classic 
Shell Theory is sufficient for fully characterizing 
wrinkling behavior on its own. 

The Mixed Interpolation of Tensorial Components 
(MITC) shell element was developed by K. J. Bathe 
and E. N. Dvorkin to combine the plate and shell 
theories in finite element analyses. 5,6,7 The stress and 
strain tensors at each element node are interpolated 
along each edge of the element, as well as through the 
element thickness. Thus, true strain energy density 
levels may be captured. This is critical, given that 
wrinkling is best described as the minimum energy 
response for accommodating in-plane strain while also 
satisfying the inextensibility principle in the transverse 
direction. 

Mathematically, wrinkles have been described in 
terms of both TFT and shell theory. For the purposes of 
this study, to be consistent with MITC element 
formulation, the principles of shell theory are used to 
characterize wrinkle states, as described in the 
following table. 

Table 1. Shell Theory Wrinkle State Criteria 

Taut an > 0 and cr 22 > 0 

Wrinkled o u > 0 and G 22 < 0 

Slack On < 0 and g 22 < 0 


Analysis 

Method 

The computational tool used to perform these 
simulations is a Finite Element Analysis Program: 
Personal Version (FEAPpv). 8,9,10 The 4-node MITC 
shell element was added by authors to the element 
library. The Newton and modified Newton methods are 
used to compute nonlinear analysis solutions. 

Model 

The 0.2m 2 membrane finite element model is 
shown in Figure 2. The mesh was discretized with 
regard to optimal aspect ratio for best performance of 
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the 4-node MITC shell element. 11 Computation time is 
of secondary importance in this study; the model 
consists of 10000 elements (0.002nT each) and 10201 
nodes. 



Fig. 2 Finite Element Model 

The membrane is constrained at 3 of 4 comers in 
ail degrees of freedom except out of piane rotation. At 
the fourth comer of the membrane, a small 
displacement is applied in the direction of the +X axis 
(5x = 0.001m) and constraints are applied against 
translation in the Y and Z directions, as well as against 
rotation about the X and Y axes. 

The magnitude of applied displacement was set 
with consideration for the assumption that this 
analytical study is based on linear elastic behavior of 
the membrane. Resultant stress levels are not to exceed 
the 5% yield point for H-type Kapton; refer to Table 2 
for assumed material properties. 


Table 2. H-type Kapton Material Properties 


Property' 

Value 

Unit 

Tensile Strength, Ultimate 

243.04 

MPa 

Yield Point (3% offset) 

95.06 

MPa 

Yield Point (5% offset) 

123.48 

MPa 

Elongation, Ultimate 

80 

% 

Young’s Modulus, E 

2959.6 

MPa 

Poisson’s Ratio, v 

0.34 


Density, p 

1.42 

g/cm 3 


Cases 

Parameters discussed in the previous section are 
held constant for all numerical simulations. Each 
analysis case performed in this study was derived based 
on the need to systematically evaluate effects of 
additional analytical variables. Table 3 lists the variable 
parameters treated. 


Table 3. Parameters Varied for Analytical Study 


Parameter 

Values 

Units 

Thickness, t 

50, 125, 150, 175, 
300, 500, 1000 

um 

Gravity' 

0 or (p-A-t)-G 

M-N 

Initial Flatness 
Imperfection 

- 1 0%t < NRN < + 10%t 

pm 


For each thickness, geometrically nonlinear 
analyses were performed in order to assess the 
wrinkling behavior both with and without the gravity 
effect, which is assumed perpendicular to the plane of 
the membrane. For the cases including this effect, 
membrane sag due to gravity is taken to be the initial 
distorted shape prior to applying the shear load (as an 
enforced displacement, 5x = 0.001m). 

For cases not including the gravity effect, the initial 
flatness of the membrane is assumed to not be perfect. 
That is, in reality, no square membrane surface can be 
considered perfectly flat when it is suspended at each of 
its 4 comers. In this study, the initial deformed shape 
(or, flatness imperfection) is determined in one of two 
ways. The first method is to apply random values of 
out-of-plane distortion at each node. These values are 
±10% of the membrane thickness and are taken from a 
normal random number (NRN) generator, and they are 
applied as initial conditions. Alternatively, the initial 
deformation is assumed to correspond to the sag shape 
due to gravity, but applied as a geometric initial 
condition, rather than taken as the result of an applied 
load. 

Table 4 summarizes each of the analysis cases 
considered in this study. All cases, but the first (“Self- 
Weight”), are considered to be wrinkling analysis cases 
for which a detailed study of the results has been 
conducted. The “Self-Weight” case is taken to be a 
baseline for comparing the quantified effects due to 
gravity simulated in the other three wrinkling cases. 
Note that each membrane thickness given in Table 3 is 
considered for all cases. 

Results 

Self- Weight Case 

The purpose of this case is to simulate the static 
effect of a (1G) gravity field on a square membrane 
held fixed (5DOF) at ail four comers. The gravity field 
is simulated by computing an effective gravity force 
(Fq) for each membrane thickness, which is evenly 
distributed over the entire area (A) and applied 
discretely at each node. 

f igure 3 illustrates the resultant deformed shape 
that is characteristic for each membrane thickness, flom 
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Z axis 


Table 4. Description of Analysis Cases 


Case 

Applied 

Load 

Initial Condition 

Boundary 

Condition 

Illustration 

“Self-Weight” 

Effective 
Gravity 
Force at 
Each Node 

None 

5DOF Constrained 
at All 4 comers 

T 

N 

♦Y 

• • 

“Gravity” 

8x=0.001m 
at one 
corner 
(X=0.2m, 
Y=0.0m) 

Restart of Self- 
Weight Case 

5DOF Constrained 
at 3/4 comers; 
4DOF constrained 
at comer of 
applied load 

SH 1 

“No Gravity 
NRN” 

5x=0.001m 
at one 
comer 
(X=0.2m, 
Y=0.0m) 

Flatness 

Imperfection due 
to random 
distribution of 
distortions (8z) 

5DOF Constrained 
at 3/4 comers; 
4DOF constrained 
at comer of 
applied load 

♦*«*« 

0 0 

“No Gravity 
S-W” 

Sag shape due to 
Self-Weight 
analysis 

T » s 

M ■}. 

♦r t 



50 pm to 1000 pm. Note that the shape of any cross 
section profile is not a catenary. 


xio* (m) 



Fig. 3 Characteristic Deformation Due to Gravity 


For comparison with the three wrinkling analysis 
cases, total strain energy density (U) for each 
membrane thickness is presented in Figure 4; a non- 
dimensional structural scale (L/t) parameter is given 
along the abscissa. 

It is clear that U rapidly decreases with thickness. 
The rate of decrease is related to the rate of decrease in 
membrane stiffness, which is a compound function of 
both t 2 and t 3 for in-plane and bending stiffness, 
respectively. Overall wurk required to distort the 
membrane is significantly less for small thickness than 
for large; therefore, this relationship is expected. 
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Fig. 4 Strain Energy Density Due to Gravity Sag 


Gravity Case 

The solution of this case is obtained from a 
nonlinear solution re-start of the “Self-Weight” case 
while applying the comer shear force. In this sense, the 
simulation may be considered to be the result of an 
initial applied load, or force. (This is contrary to the 
“No Gravity S-W” case, described below, in which the 

giaviiy cnoci is ao a gcumtuw niiucu 

condition, or distortion.) 

Although simulations for each thickness given in 
Table 3 w'ere performed, for brevity only the deformed 
shapes for t = 50, 175, 300, and 1000 pm are illustrated 
in Figure 5. Simulation results indicate that when t < 
500pm, maximum out-of-plane distortions (5z, in either 
the +Z or -Z direction) decrease as thickness decreases. 
However, there is no clear trend for maximum stress 
levels that can be related to thickness; the difference 
between maximum and minimum values is only 5.7%. 

Element-by-element analysis of wrinkle state, as 
defined in Table 1, for each membrane thickness has 
revealed that for t > 300pm, wrinkle criteria based on 
major principal stresses does not accurately characterize 
membrane distortion because some regions do carry 
some level of compressive stress without buckling due 
to sufficiently high levels of bending stiffness. This 
may be confirmed by examining the deformation modes 
given in Figure 5; true wrinkling behavior (a result of 
localized buckling) is dominant only when t < 300pm. 

No Gravity NRN Case 

This analysis case was prepared by enforcing, as a 
geometric initial condition, small out-of-plane (5z) 
distortions at each node location. The magnitude of 
these distortions is not greater than 10% of the 
membrane thickness, and they were taken from a 
normal random number (NRN) generator. 


t = 50 pm 
Max. |5z| = 2.6 mm 
CTpmax = 141.4 MPa 
Max. |Xxy| = 55.8 MPa 

*'1 

•v.axtmi 

♦X ax* (m) 

t= 175 pm 
Max. |Sz| = 2.9 mm 
opmax = 142.8 MPa 
Max. |Xxy| = 56.9 MPa 

1 m *Y »xi» (m) 

♦x m frn) *** 

t = 300 pm 
Max. |5z| = 3.0 mm 
CT Pmax = 144.5 MPa 
Max. |x xy | = 57.5 MPa 

: - - jj 

“ ■; ” 

* ^ -y in (m) 

♦X«xw(m) 

t = 1000 pm 
Max. |5z| = 2.6 mm 
CTpmax = 150.0 MPa 
Max. IXxyl = 59.1 MPa 

- ’1 s ’ 

• y" ,■ *Y txu (m) 

/'•* 

•X *** (m) V, 

9> • 


Fig. 5 Wrinkled Membranes Due to Applied Corner 
Shear Load, Including Gravity Field 


This initial condition was chosen in order to induce 
a perturbation in the flatness of the membrane - in the 
absence of a gravity effect. The perturbation is 
necessary because the assumption of a perfectly flat 
surface is an inherently stable system under static load 
conditions and, therefore, cannot be numerically 
induced to wrinkle. 

Other analytical methods have enforced, as an 
initial condition for wrinkling simulations, the primary 
eigenmode of deformation for the given structure. 
How'ever, this can arguably bias the solution in favor of 
deformation modes that correspond only to trajectories, 
or eigenvectors, related to that particular eigenmode. 

As in Figure 5 for the Gravity case, Figure 6 
illustrates wrinkled membranes only for t = 50, 175, 
300, and 1000 pm. A quantitative comparison of the 
stresses and deformations associated with these results 
and those given in Figure 5 reveals no significant 
differences. However, there are obvious qualitative 
differences as thickness decreases. 
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t = 50 pm 
Max. |8z| = 2.4 mm 
tfPmax= 141.1 MPa 
Max. |x xy | = 55.7 MPa 


*X Brt« (m) f -*» 

t = 175 pm 
Max. |8z| = 2.8 mm 
Cpmax = 142.8 MPa 
Max. |Xxy| = 56.9 MPa 


(mj 

♦X*x1«(m) *■** 

t = 300 pm 
Max. |8z| = 2.8 mm 
Clpmax = 144.7 MPa 
Max. |x xy | = 57.6 MPa 

* 

»X axis (m) "1 

t = 1000 pm 
Max. |Sz[ = 2.6 mm 
tfpmax = 143.8 MPa 
Max. x xy | = 57.8 MPa 

i 

f:. 

♦x txi* (m) *’* 


Fig. 6 Wrinkled Membranes Due to Applied Corner 
Shear Load, Including Gravity Field and Initial 
NRN Imperfection 

Therefore, in order to truly characterize wrinkling 
behavior it is necessary to quantify wrinkle orientation 
in addition to amplitude. This has been done, for this 
study, in a rudimentary way by plotting maximum and 
minimum principle stresses at a constant cross section 
in the membrane, for each thickness^ as shown in 
Figure 7. The most significant detail to notice from this 
comparison is that minor principle stresses tend to 
reverse, and with greater frequency, as thickness 
decreases. This is the cause of increasingly numerous 
and skewed wrinkles when thickness is less than 300 
pm. 

No Gravity S-W Case 

This wrinkling case differs from those presented in 
the previous two sections in that the assumed initial 
condition is determined by the appropriate degree of 
gravity sag deformation, depending on thickness, that 
was computed in the Self-Weight case. This is unlike 
the Gravity analysis case because the initial condition is 


1 This topic is expanded and thoroughly treated in a 
separate study entitled “A Spectrum Evaluation Method 
for Wrinkled Membranes Using Two-Dimensional Fast 
Fourier Transform” by T. Iwasa, ct. al. 3 



Fig- 7 CT n and ct 22 Along Membrane Cross Section: 
Y = 0.14m 


applied as a geometric constraint, rather than as a 
distributed applied load (F 0 ). 

Figure 8 gives analysis results similar to those 
reported in Figures 5 and 6, for t = 50, 175, 300, and 
1000 pm. These results reveal similar trends with 
respect to wrinkle state; when t > 300 pm deformation 
may not be accurately characterized by the term 
“wrinkled”, according to the state definitions given in 
Table 1. In addition, stress levels and out-of-plane 
deformation are not significantly different than the 
other two wrinkling cases. 

However, the deformed state when t = 50 pm is 
different than the previous two cases. In addition, 
differences are discovered when carefully analyzing 
membrane deformation at constant cross sections using 
the wrinkle scale parameter (W/t). This is discussed in 
greater detail in the following section. 

Discussion 

In order to visualize how wrinkling behavior is 
effected by the different initial conditions it is 
beneficial to examine cross section plots at constant 
thickness for each analysis case. The non-dimensional 
wrinkle scale parameter W/t is plotted versus structural 
scale L/t, where W=amplitude and L=edge length. For 
illustration, here, the cross section plots for t = 50, 175, 
300, and 1000 pm at y = 0.14m arc shown in Figure 9. 
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CTPmax = 138.5 MPa 
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• « 



Fig. 8 Wrinkled Membranes Due to Applied Corner 
Shear Load, Including Gravity Field and 
Self-Weight Gravity Sag 


A detailed examination of such plots for each 
membrane thickness reveals that deformation mode is 
independent of initial conditions when t < 150pm. Also, 
for any thickness, the deformation modes for “Gravity” 
and “No Gravity S-W” are the same, but offset by the 
initial 8z offset due to the gravity sag shape. The 
disparity between these two deformation modes and 
that for the “No Gravity NRN” case becomes noticeable 
when 175 pm < t < 1000pm and a phase shift between 
wrinkle peaks becomes apparent. But, when t = 
1000pm, deformation mode is again independent of 
initial conditions. Another prominent trend to notice is 
that wavelength increases with thickness, regardless of 
assumed initial conditions. 

Figure 10 illustrates the relationship between the 
wrinkle (W/t) and structural (L/t) scale parameters. It is 
important to note that this relationship is linear for the 
Gravity case. This indicates that analysis results may be 
scaled for correlation to ground-based experimental 
data - but only for square geometries. That is, previous 
research studying similar behavior but for a circular 
membrane has revealed that the same relationship is not 
linear. 12 Rather, it increases nonlinearly towards a 
constant value as L/t increases. So, the relationship 
between wrinkle scale and Structural scale depends on 
geometry/, not initial conditions. 



Fig. 9 (W/t) Along Membrane Cross Section: 
Y = 0.14m 


Overall Maximum W/t vs L/t 
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Fig. lft Overall Maximum (W/t) versus (L/t) 
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Plotting total strain energy density versus structural 
scale reveals tire same relationship for all three 
wrinkling analysis cases as was revealed for the Self- 
Weight case. Refer to Figure 11, which shows that 
volumetric strain energy decreases rapidly with 
thickness and approaches a constant value. Therefore, 
the minimum energy response for accommodating in- 
plane strain (in the form of wrinkles) approaches a 
constant as thickness becomes small. However, as 
previously noted, when t > 500pm wrinkling is not an 
accurate description of the deformation due to in-plane 
shear. 



Fig. 11 Strain Energy Density versus (L/t) for All 
Wrinkling Cases 


With respect to strain energy density, it is also 
interesting to note that the ratio of results for the 
Gravity and Self-Weight cases increases rapidly and 
approaches a constant value as thickness decreases 
(refer to Figure 12). 


Strain Energy Density Ratios: Uorivtcy/Usw 
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Fig. 12 Strain Energy Density Ratio: 

^Gravity/ Uself-Weight 

Additional analyses have revealed that the ratio of 
strain energy density in the Gravity case to that for 
either of the other two wrinkling cases is 1.0 ±5%, 
regardless of thickness. Therefore, initial conditions do 
not significantly impact (lie energy slate of a wrinkled 
membrane. Other parameters relating wrinkle amplitude 
to orientation are necessary to hilly quantify and qualify 
wrinkle behavior. 


Summary 

A systematic analytical study of the affects of 
thickness, gravity, and initial conditions on membrane 
wrinkling behavior has been presented. The universally 
applied static load in this study is in-plane shear, 
applied at one comer of a square membrane held fixed 
at the remaining three comers. This study is unlike 
others previously conducted on wrinkling behavior 
because of the geometry and numerical methods used to 
simulate behavior. 

Detailed analysis of the simulation results for three 
different wrinkling cases reveal important relationships 
between deformation and thickness, as well as initial 
conditions: 

• As thickness decreases, wrinkle behavior is 
increasingly sensitive to assumed initial 
conditions. 

• As thickness increases, greater than 300 pm, 
membrane deformation may no longer be 
mathematically characterized as wrinkling. 

• Overall stress and strain levels are not 
significantly effected by assumed initial 
conditions, though wrinkle orientation is. 

• Though the magnitude of out-of-plane 
distortion is nearly constant regardless of 
thickness or initial condition, wrinkle 
wavelength increases with thickness for any 
assumed initial condition. 

• By plotting wrinkle scale with respect to 
structural scale, it is shown that wrinkle scale is 
independent of the nonlinear relationship 
between thickness and stiffness. 

• The relationship between wrinkle scale and 
structural scale is linear for the Gravity case in 
this study, though it is nonlinear for a circular 
membrane subjected to in-plane torsion. 
Therefore, this relationship depends on 
geometry, not initial conditions. 

• The ratio of overall strain energy density in the 
Gravity case to that for either of the other two 
wrinkling cases is approximately 1.0. 

All of these results may be useful for correlating 
test data from a ground-based experiment. In addition, 
they reveal the fact that quantifying wrinkle behavior 
due to in-plane shear is a laborious and time-consuming 
task when utilizing finite element simulations. Because 
knowledge of precise wrinkle characteristics such as 
orientation and wavelength may be necessary for 
optical applications of membrane structures, this 
analytical study may be useful, as the background, for 
developing alternative wrinkling analysis tools capable 
of predicting such details outright. 
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